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Splitting and oscillation of Majorana zero modes in the p-wave BCS-BEC evolution 

with plural vortices 
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We investigate how the vortex-vortex separation changes Majorana zero modes in the vicinity of 
the BCS-BEC (Bose-Einstein condensation) topological phase transition of p-wave resonant Fermi 
gases. By analytically and numerically solving the Bogoliubov-de Gennes equation for spinless p- 
wave superfluids with plural vortices, it is demonstrated that the quasiparticle tunneling between 
neighboring vortices gives rise to the quantum oscillation of the low-lying spectra on the scale of the 
Fermi wavelength in addition to the exponential splitting. This rapid oscillation, which appears in 
the weak coupling regime as a consequence of quantum oscillations of quasiparticle wave functions, 
disappears in the vicinity of the BCS-BEC topological phase transition. This is understandable from 
that the wave function of the Majorana zero modes is described by the modified Bessel function 
in the strong coupling regime and thus it becomes spread over the vortex core region. Due to the 
exponential divergence of the modified Bessel function, the concrete realization of the Majorana 
zero modes near the topological phase transition requires the neighboring vortices to be separated 
beyond the length scale defined by the coherence length and the dimensionless coupling constant. 
All these behaviors are also confirmed by carrying out the full numerical diagonalization of the non- 
local Bogoliubov-de Gennes equation in a two dimensional geometry. Furthermore, this argument 
is expanded into the case of three-vortex systems, where a pair of core-bound and edge-bound 
Majorana states survive at zero energy state regardless of the vortex separation. 

PACS numbers: 05.30.Fk, 03.75.Hh, 03.75.Ss, 74.20.Rp 



I. INTRODUCTION 

A Majorana fermion is a relativistic particle equiva- 
lent to its anti-particle, which was originally proposed by 
Ettore Majorana in 1937 [l|. It has been believed that 
the Majorana zero modes offer the key to an understand- 
ing of the neutrino mass problem [2J and a fundamental 
building block of topological quantum computing [3|t51 ■ 
Recently, it has been, moreover, predicted that they are 
hidden in superconducting and superfruid materials [3j. 

A spinless chiral p-wave supernuid may be an effec- 
tive model that describes the low energy properties of 
various p-wave superfluids without time reversal symme- 
try, such as p-wave resonant Fermi gases 043, half- 
quantum vortices [l7H22j |. noncentrosymmetric super- 
conductors |23rf26| , and supernuid-ferromagnet insulator 
junctions formed on the topological insulator [27j. Using 
this model with singular vortices, the low energy excita- 
tions consist of two characteristic quasiparticles, such as 
the vortex core bound state and edge bound state. Their 
eigenenergies in a chiral k x — ik y state with arbitrary 
vortex winding number k are found to be prop ortional 
to an azimuthal quantum number !eZ |, Qj], |28| - |33 |, 

Ej, c '° ' (x£ — ^^i, where the weak coupling BCS regime is 
assumed. The noticeable consequence is that the lowest 
eigenenergy of E^ c can be exactly zero when k is odd. 
The zero energy quasiparticle is composed of the equiv- 
alent contributions from the particle and hole, and thus 
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its creation is describable with a sclf-Hcrmitian operator 
if E _ = r)E=o, called the Majorana zero modes. 

The remarkable fact arising from the self-Hermitian 
property is that the Majorana zero modes and their 
host vortices obey neither Fermi nor Bose statistics, be- 
cause they break the ordinary anti-commutation relation 
nj^ = 1/2 and { Ve ^, Ve ^} = 1 0, Q3, M, M, M, H ■ 
An ordinary fermion can be restored by taking account 
of linear combination of two Majorana operators, called 
the complex fermion. The p-wave supernuid with well- 
separated 2n vortices may contain 2n-th Majorana states, 
leading to the 2™-th degeneracy of the many-body ground 
state differentiated by the occupation number of the zero 
energy complex fermions. Hence, the degenerate quan- 
tum state can be a topologically protected qubit [33, [38| , 
whose manipulation can be implemented by braiding vor- 
tices as a consequence of the non-abelian statistics of vor- 
tices. For instance, a discrete set of the unitary group 
which manipulates the qubit can be im plem ented by the 
braiding operation of four vortices [171 |39j . In addition, 
it has recently been proposed that the continuous_ manip- 
ulation can be realized in three- vortex systems [401 ] . 

While this system may offer the promising method of 
the fault-tolerant quantum computation [3 5], it has re- 
cently been revealed that the intervortex tunneling and 
thermal fluctuations of vortic es g ive rise to the decoher- 
ence of the topological qubit [4l| . Since the zero energy 
wave function is localized within the core radius when the 
vortices are well separated from each other, the quasipar- 
ticle tunneling between the Majorana zero modes lifts 
the degeneracy of the many-body ground states expo- 
nentially with respect to the ratio of the vortex distance 



and the core radius [41lT43j . This gives rise to the de- 
cohercncc of the Majorana qubit and may be critical for 
the implementation of a fault-tolerant operation. This 
exponential splitting has been found to be ubiquitous 
in various systems associated with non-abelian anyons, 
such as the non-abelian quasiholes of the v= | fractional 
quantum Hall state [44.|45|. Kitaev's honeycomb lattice 
model [46], and the generic anyon model [43|. In ad- 
dition to the splitting of Majorana zero modes, for the 
weak coupling limit, the rapid oscillation of the eigenen- 
ergy due to the quantum nature of the Majorana zero 
mode appears on the scale of the Fermi wavelength with 
varying the vortex separation |4l|. Here, we expand the 
argument into the vicinity of the BCS-BEC topological 
phase transition (TPT) beyond the weak coupling limit, 
which can be driven by a magnetic Feshbach resonance 
in p-wave resonant Fermi gases [a-KUfl, [2a . l48l - [53j . 

The aim in this work is to clarify an unsettled ques- 
tion, that is, how the tunneling between the neighboring 
vortices changes the Majorana zero modes in the BCS- 
BEC evolution of p-wave resonant Fermi gases. In the 
limit that neighboring vortices are sufficiently separated 
from each other, the Majorana zero mode for an odd vor- 
ticity vortex survives until the BCS-BEC TPT point is 
approached from the weak coupling BCS limit. In the 
BEC regime beyond the TPT point, the low- lying quasi- 
particle spectrum becomes trivial. Recently, Gurarie and 
Radzihovsky [32[ have revealed that in contrast to the 
BCS limit, the wave function of the zero energy states 
in the vicinity of the TPT is described by the modified 
Bessel function, leading to the derealization beyond the 
vortex core region. This has been confirmed by our pre- 
vious study based on the self-consistent calculation [13j |. 

In this article, we examine the splitting and oscillation 
of the Majorana zero modes through the inter- vortex tun- 
neling in the vicinity of the TPT, based on the analyti- 
cal and fully numerical calculations of the Bogoliubov-de 
Gennes (BdG) equation. It is demonstrated that in the 
weak coupling limit, the splitting of Majorana zero modes 
due to the inter-vortex tunneling is characterized by a 
single dimensionless parameter composed of the ratio of 
the coherence length and vortex separation. It turns out 
to depend on the additional length scale defined by the 
coherence length and the dimensionless coupling constant 
in the vicinity of the TPT. Due to the exponential diver- 
gence of the zero energy wave function, the realization 
of the Majorana qubit requires the neighboring vortices 
to be separated from each other beyond the new length 
scale. Furthermore, a new topological qubit recently pro- 
posed by Ohmi and Nakahara 40] can be continuously 
manipulated by braiding of three vortices, in contrast 
to a discrete set of the unitary group in four- vortex sys- 
tems [l3,|39|. Nevertheless, the stability of the Majorana 
zero modes in three- vortex systems has never been stud- 
ied so far. Hence, we expand the argument in two-vortex 
systems into three- vortex systems, where a pair of edge- 
and core-localized Majorana modes is found to always 
survive regardless of the vortex separation. 



This paper is organized as follows. In Sec. II, we de- 
scribe our theoretical framework based on the BdG equa- 
tion with a non-local potential in two dimensional ge- 
ometry. Here, since the BdG equation reduces to an 
eigenvalue problem with a huge matrix, we shall out- 
line here about the numerical diagonalization with the 
discrete variable representation and Krylov subspace it- 
eration. This is also supplemented in Appendices A and 
B. After the Majorana solution of the BdG equation is re- 
viewed in Sec. Ill A, the expression on the splitting of the 
Majorana zero modes due to the inter- vortex tunneling is 
analytically derived and compared with the full numer- 
ical calculation. The details on the splitting modes are 
presented and discussed throughout the remaining part of 
Sec. III. In Sec. IV, we expand our study into three- vortex 
systems. The final section is devoted to conclusions and 
discussion. In addition, we give in Appendix A supple- 
mentary information concerning the issue on the com- 
plex eigenvalues which appear when the BdG equation 
within low-energy approximation is numerically solved. 
In Appendix B, we describe the details on the numeri- 
cal diagonalization of the non-local BdG equation in two 
dimensional geometry based on the discrete variable rep- 
resentation I54rl56i. 



II. THEORETICAL FORMULATION 

A. Bogoliubov-de Gennes equation and p-wave pair 
potential 

Here, we consider spinless fermions interacting via an 
effective pairing potential V(r\, r 2 ) with the mass M. It 
is convenient to introduce a spinor in the Nambu space, 
^(ri) = [^/'(ri),^'''(ri)] T , with the creation and annihila- 
tion operators of spinless fermions i/j' and i\). Using the 
definition of the pair potential, 



A(r 1 ,r 2 ) = ~V(r 1 ,r 2 )(ip(r 1 )^(r 2 )), 



(1) 



the Hamiltonian within the mean-field approximation is 
given by 



H = E Q + - J dn /dr 2 * + (ri)/C(ri,r 2 )*(r2) 



where the matrix K, is given as 



/C(n,r 2 ) = 



H ( ri )S( ri - r 2 ) 
-A*(n,r 2 ) 



A(n,r 2 ) 

-H*( ri )6( ri - r 2 ) 



(2) 



,(3) 



with Hq(v) = —-^jj — M- Throughout this paper, we set 



% = ks = 1. The p-wave pair potential must satisfy the 
symmetry on the orbital degrees of freedom, 



A(r 2 ,n) = -A(ri,r 2 ). 



(4) 



The mean-field Hamiltonian in Eq. ([2]) can be diago- 
nalized by introducing the unitary transformation to the 



quasiparticle basis with r\ v = [77,,, ?/J] T , 
*( r i) =^2'G>u{r 1 )ri u , 



(5) 



where u v is a 2 x 2 matrix and the matrix elements de- 
scribe the quasiparticle wave functions. This is required 
to satisfy the orthonormal condition, 



ui(ri)u V '(ri)dri = S v , 



v' 1 



(6) 



and completeness conditions, ^„ u y (ri)uj,(r 2 ) = S(r\ — 
r 2 ). Also, the fermion operators r\ v and r\ J obey the anti- 
commutation relation, {t?i/)%/} = $i/,i/' and {?7iy,?/i/'} = 
{^J,?7^,} = 0. The mean-field Hamiltonian in Eq. ([2} is 
now diagonalized in terms of this basis with the quasi- 
particle energy E v as W = E + \ J2 E v rjlt3rj v . This di- 
agonalization leads to the Bogoliubov-de Gennes (BdG) 
equation, j dr%K.(r\, r 2 )u„(r 2 ) = E v u v (r{). Here, it is 
found that the matrix elements of u yield (u)22 = (w) l:L 
and (u)i2 = (11)21, because of the symmetry, 



-fi/C*(ri,r 2 )fi =/C(ri,r 2 ) 



(7) 



where ti j2j 3 denote the Pauli matrices. Hence, the BdG 
equation reduces to the equation for the quasiparticle 
wave function u u = (u)n and v v = (u)2\, 



dr 2 K(ri,r 2 ) 



Uy(r 2 ) 
v v (r 2 ) 



E u 



u v (r x ) 
v v {r x ) 



(8) 



This equation describes the energy eigenstate in the pres- 
ence of the non-local pair potential A(ri, r 2 ). 

In order to clarify low-lying quasiparticle structures 
in chiral p-wave superfluids, we directly solve the BdG 
equation ([5} with the non-local pair potential A(r 1 ,r 2 ). 
First, we derive the general expression of A(ri,r 2 ) in 
a three-dimensional coordinate r = (x,y,z). The pair 
potential A(ri, r 2 ) is expanded to the Fourier series with 
respect to the relative coordinate r\ 2 = v\ — r* 2 , 



n f dk . 

A(r - r2 W W A 



r\ +r 2 



k 



Ah-r\ 



(9) 



with the relative wave vector k = (k x ,k y ,k z ). Here, we 
assume A(r, k) to be expanded as 



A(r,fc) = J2 A m (r)T m (k). 



(10) 



m=0,±l 



The function A m (r) on the center-of-mass coordinate 
r describes the spatial variation of the pair potential 
around vortex cores, as we will see in Sec. II B. The func- 
tion r m (fc) in Eq. (TlOl) is a symmetry factor that de- 
scribes an attractive interaction in p-wave channel. In 
the uniform system, the symmetry factor is obtained by 
replacing the interparticle potential V in Eq. ([T]) to a 
model potential HEllO] 



I m\l&) — 



with 
2 + k 2 o 



(11) 



Here, k m is the eigenstates of the angular momentum 
operator of L = 1 in relative coordinate, L z k m = mk m 
and k±i = ^f(k x ± ik y ) and kg = k z . The parameter ko 
in Eq. (| 1 1 j) corresponds to the inverse of an effective in- 
teraction range. Since the diluteness of systems requires 
fc^ 1 -C fcp 1 , the low energy physics of the BdG equa- 
tion (JSJl under A(r,k) within k~k~F may be describablc 
with T(k)fHj-k m . Note that the symmetry factor yields 
r m (|fc| — >oo) =0 in the high energy limit. For the numer- 
ical diagonalization of the BdG equation (J8J), however, we 
further replace the expression of T m (k) in Eq. (|lll) to 



r m (fe) 



k 
k 



-(k 2 -kl)/k 2 i 



(12) 



Here, the additional parameter k\ is introduced, which is 
assumed to be an order of the Fermi wavelength ki — kp. 
Note that this expression on T m (k) in Eq. (|lll) correctly 
reproduces the low- and high-energy behaviors of the 
original form in Eq. dill) . As we will see below, in particu- 
lar, the behavior r m (|fc| — >oo) =0 is crucial for preserving 
the p-wetve symmetry of A(t*i , r 2 ) in Eq. (U), which guar- 
antees the eigenvalues of Eq. (|SJ) to be real in numerical 
diagonalization. 

By substituting Eqs. ([ID]) and (QJ) into Eq. © and per- 
forming the integral over k, the pair potential A(ri,r 2 ) 
in a real coordinate is given by 



A(ri,r 2 ) 



E 

m=0,±l 



A; 



ri + r 2 



F m (ri 2 ) (13a) 



r±i(r 12 ) = T 



r (n 2 ) 



^0 
16tt 3 / 2 



(x 12 ±iyi 2 )e 



Ingi** 



-;.- 



16tt 3 / 2 



^i 2 e 



(13b) 



(13c) 



It is easy to see that the eigenenergy of the BdG equation 

© is given by E(k) = y/e 2 k + \ YZ T m {k)A m \ 2 with e k = 
,2 
-ryjj — /i, when the uniformity of the pair potential A m is 

assumed. 

Throughout this paper, we assume the chiral p-wave 

pairing state, (A+i,Aq,A-i) = (0,0, .4). This can be 

realized in a two-dimensional geometry, where the pair 

potential and wave functions reduces to A(r) —A(x,y), 

u v (r) — u u (x,y) and v v (r) — v v (x,y). Then, the BdG 

equation ((SJ in the two-dimensional geometry is rewritten 



dp 2 K,(pi,p 2 ) 



U V (P2) 
V V (P2) 



= E V 



u v (pi) 

Vu(Pl) 



(14) 



where p = (x,y). The pair potential A(pi,p 2 ) = 
J dz\ J dz 2 A(r\, r 2 ) in a two-dimensional plane is given 
from Eq. ([13]) by 



A , n , ixi2 + yi2 - lP12 i 2fc ° +|# , ( p\ - p-i 



87rfc ' 



,(15) 



where p\2 = (£12,2/12) = Pi — Pi- It is obvious that the 
expression of the non-local pair potential in Eq. (|13[) and 
(fT5f still satisfy the p-wave symmetry in Eq. (j4]) and the 
Hamiltonian density /t(r 1) r 2 ) is Hermitian, which guar- 
antees the eigenvalue E„ to be real. 

Since this article focuses on the low-energy quasipar- 
ticles, it might be convenient to carry out the approx- 
imation with fco — > 00 hi Eq. (1121) . which reduces the 
symmetry factor to r m (fc) w T-k m . This simplifica- 
tion, which neglects the effective size of the Cooper pair, 
changes the the non-local pair potential to the local ex- 
pression A(ri,r 2 ) « S(n - r 2 )^J2m^),±i{ A ( r )i p m } 
in derived in Eq. (|A1|I . Here, P m consists of the linear 
combination of spatial derivatives as described in Ap- 
pendix A. It is important to note that the A(ri,T"2) re- 
sulting from fco — > 00 no longer satisfies the p-wave sym- 
metry in Eq. Q. Then, the non-local BdG equation 
© reduces to the local form described in Eq. (JA3I) with 
fc{ri,r 2 ) = K.(r)S(ri — r 2 ), where the matrix K.(r) re- 
sults in the non-Hermitian matrix. As we will discuss in 
Appendix A, this non-Hermitian matrix may contain the 
complex eigenvalues, especially for E v — > 0. In practice, 
a non- vanishing imaginary part of the eigenvalues E v , 
which appears when the non-Hermitian matrix is numer- 
ically diagonalized, gives rise to the abrupt jump of the 
corresponding real part. 

It is known that non-Hermitian matrices also appear 
in the BdG equation in Bose-Einstein condensates, which 
describes a small fluctuation around a given ordered state 
within the mean-field approximation for dilute Bose sys- 
tems [53, [5g|. In the context of Bose-condensed systems, 
the physical meaning of the complex eigenvalues has been 
discussed by a large number of authors [59I462J and has 
been found to be associated with the dynamical instabil- 
ity of a given ordered state without any dissipation. In 
contrast to Bose systems, however, the non-Hermitian of 
1C(r) is in consequence of the approximation on A(t*i, r 2 ) 
at fco — >oo, which no longer satisfies the p-wave symmetry 
in Eq. (U), and the imaginary value of the eigenenergies 
must be unphysical. Hence, the BdG equation (|14p with 
non-local pair potential A(ri, r 2 ) in Eq. (|15[) has to be di- 
rectly solved in order to avoid the emergence of complex 
eigenvalues in numerical diagonalization. The non-local 
BdG equation (|14[) can be numerically solved by using 
the discrete variable representation (DVR) and Krylov 
subspace iterative method as described in Appendix B. 



B. Vortex configuration and BCS-BEC evolution 

Throughout this work, we assume the spatial shape of 
the pair potential around vortices expressed through A 



Here, Kj is the winding number of the j-th vortex, 
Pj = p — Rj denotes the coordinate centered at the j-th 
vortex core Rj, and 9j = ta.n~ 1 (yj/xj) is the polar angle. 
The vortex positions {Rj}j=i,— ,n w are set to be i?i = 
D v (^,0) and R 2 = D v (— 5,0) for N v — 2 systems and 

R 1=J D v (-^,0), R 2 =D v (-±,&), i?3 = £> v (-I,-^) 

for N v — 3 systems (see also Fig. [1]). The BdG equa- 
tion (fl4|) with Eqs. (fT5|) and (fT6|) gives the low-energy 
quasiparticle spectra for spinless p-wave superfluids with 
plural vortices, where the BCS-BEC evolution is parame- 
terized by varying the strength of Ao and fx, as described 
below. Here, we impose the rigid boundary condition on 
the quasiparticle wave functions at the radius R= 150fcp 

as u n {\p\=R)=v n (\p\=R)=0- 

The resulting BdG equation © contains four length 
scales, such as the mean interparticle distance fcp , the 
coherence length £ = (2-Ep/Ao)fcp , the distance between 
neighboring vortices D v = \Rj — i? J+ i|, the radius of 
the system R. Here, since we are interested in the sit- 
uation, fc^ 1 < fcp 1 < £ <C D v < R, the large number 
of the DVR basis N has to be taken. In cold atoms 
confined in a harmonic trap, for instance, D v /£ is esti- 
mated as D v /£ < Rtf/Z ~ VN at om/k F Z ~ 0(100) with 
the Thomas-Fermi radius of the cloud i?TF and the to- 
tal particle number V atom = O(10 5 ). In this calcula- 
tion, we set the number of the DVR basis function to 
be N = 300, which requires the large size of the mem- 
ory. To overcome this issue, the huge and dense matrix 
is numerically diagonalized with the shift-invert Lanczos 
algorithm [63] . This algorithm can reduce the eigenvalue 
problem to the 2JV 2 -dimensional linear equation, which 
is iteratively solved by the Krylov subspace method, such 
as the generalized minimal residual method with a pre- 
conditioner 64]. 

It is convenient to introduce the dimensionless param- 
eter 



fc M £o 






ku 



fcp 



m ^F A 

fc 



(17) 



which can parametrize the BCS-BEC evolution in p-wave 
resonant Fermi gases [13| . Here, fc M = y/2M\/j,\ is defined 
with the chemical potential fi. The length scale fc^ 1 is 
found to describe the shortest wavelength of the low- lying 
quasiparticle wave function. In the weak coupling limit 




A' v 



A{p) = A \{ e iK ' Si tanh 



j=i 



£0 



(16) 



FIG. 1: Schematic picture about the vortex configuration in 
the case of N v = 2 (a) and N v = 3 (b) with the core radius £. 



where p ss Ep 3> Ao, it turns out to be equivalent to 
kfj,£,o ~ &f£ 3> 1- As the amplitude of the pairing inter- 
action increases, however, the values of p and A deviate 
from those in the weak coupling regime. In practice, p 
(A) decreases (increases) gradually with increasing the 
pairing interaction and becomes zero at certain value of 
the coupling constant [g, 0, |9h111 [l3[ . The point at which 
p becomes zero (k li £,o = Q) is regarded as the TPT point, 
which is known to give rise to the drastic change of the 
low-lying spectra of chiral p-wave supernuids. The low- 
lying spectra in the regime with p>0 (fc„£o >0) is sensi- 
tively affected by the non-trivial topological structure of 
the pair potential, while the BEC regime (p < 0) beyond 
the TPT point has quasiparticle excitations with a trivial 
gap sw |/i|. Hence, we focuses on the regime with p > 0, 
i.e., fc M £o>0, in this article. 



III. SPLITTING OF MAJORANA ZERO MODES 
IN TWO- VORTEX SYSTEMS 

A. Majorana zero modes 

The odd parity of the pair potential, A(ri,r 2 ) = 
— A(r2,r , i), leads to the symmetry of the BdG matrix 
^(ri,r 2 ) described in Eq. ([7]). This gives rise to fact 
that if [ue,ve] t is an eigenstate of the BdG equation 
(|8|) with a positive energy E > 0, the corresponding 
[v,-e,v~e] T = [ v e' u e] T must be an negative energ 



the relativistic Dirac Hamiltonian[65], |66[ and the qua- 
siclassical analysis of the p-wave BdG equation. [34[ The 
zero energy solution in the BdG equation (fi"4")) is derived 
from the reduced BdG equation (IA3|) . corresponding to 
the low energy approximation of the BdG equation f| 14[) 



with r r „(fc) f=a v-k m . Following the procedure made by 



Gurarie and Radzihovsky [31 
is given as 



the zero energy solution 



ipe(p) = Afe i 



ff 



flip) 

L(p)e 



-i(k-I) 



-p/io 



(20) 



where the azimuthal quantum number I € Z for the zero 
energy state is determined by the vorticity k of the pair 
potential A-\(p) = A e iKe and A+\ = Aq = 0, through 
£=(k— 1)/2€Z. The function ft is expressed as 



flip) = h pk M l- 



1 



(fc^o) 2 



for fcu£o > 1 and 



flip) = h pk^ 



(fc^o) 



(21a) 



(21b) 



for k^o < 1. Here the parameter fc M £o is introduced 
in Eq. (fT7|) . In addition, Ji and -Z^(z) are the £-th or- 
der Bessel function and modified Bessel function and 
Af is the normalization constant. For the BCS limit 

eigenstate^ Using this Symmetry with Eq. ©, the quasi- with k ^o ~ k F £ > 1, the zero energy wave function in 

Eq. (|21a|) consists of the quantum oscillation on the scale 
of the Fermi wavelength due to Ji(pkp) and the expo- 
nential decay factor as ui~e~ p ^° cos(pkp). In contrast, 
it is found from Eq. (|21bj) for fc^o < 1 that the quantum 
oscillation of the wave function disappears and the factor 
e -p/?o m Eq (|21bj) is canceled out by the modified Bessel 



particle creation operator ry E is then found to be equiva- 
lent to the annihilation of quasiparticle with the negative 
energy, i.e., 



V-e 

'l-E 



V-E 



* 



U-E 



* T y = 



«E 

'IE 



,(18) 

L ' — -C/ J V L ~ J L ~ J J L »-" J 

where ue 1 ^ = J u E {r)^>(r)dr and v E ^ = ^ v E {r)^{r)dr. 
This consequence implies that if the quasiparticle energy 
is exactly zero, the eigenfunction always satisfies 



function for a large p as f^o(p) 
A 



VP l 



-0—X)p/£o where 



n 



us=oir) 
VE=oir) 



UE=o(r) 
vs=oir) 



(19a) 



and its operator yields the self-Hermitian condition 



V E =o =Ve=o- 



(19b) 



Since the positive and negative energy states appear as a 
pair and the total number of all eigenstates in Eq. © is 
also even, the zero energy state with u n (r) =w*(r) and 
E n = must appear as a pair [321 ] . 

Let us now revisit the zero energy solution of the BdG 
equation. It is demonstrated that there exists only a sin- 
gle zero energy Majorana state for an odd vorticity k 
vortex and none for an even vorticity, as has been re- 
vealed in the BCS limit |31j . the more generic situation 
|32j |. and the full numerical calculation [I3]. This is con- 
trast to the index theorem for zero energy eigenstates of 



■Jl — (fc^Co) 2 - This implies that the wave function is 
extended over the coherence length, which manifests the 
topological phase transition at fc^o - ^0. 



B. Splitting of Majorana zero modes 

Figure [5] is one of our main results, where we summa- 
rize the splitting of the lowest eigenenergies as a function 
of the vortex separation D v /^q. The symbols represent 
the full numerical calculation of the BdG equation (TTJ} 
with the non-local pair potential for various values of 
^/jfo, such as k^o = 10, 6.67 and 4 for the weak coupling 
regime and fc^£o = 0.63, 0.77 and 0.89 for the strong cou- 
pling regime. The solid and dashed curves depict the 
splitting energy analytically derived with Eq. (|2"Tj) . as we 
will describe below in details, e.g., in Eqs. (|2"5)) and (j2"9"j) . 
It is found from Fig. [2] that the Z? v -dependence of the 
lowest eigenenergies for the weak coupling regime within 
k^o > 1 is independent of the values of £; M £o- The expo- 
nential splitting in this regime is expressed by the single 



dimensionless parameter D v /£o, where £o reflects the lo- 
calized area of the Majorana zero mode at each vortex 
core. It will also be demonstrated analytically and nu- 
merically in Eq. ([25)1 and in Sec. Ill C that the lowest 
eigenenergies yield the rapid oscillation cos(/cf-Dv) in ad- 
dition to the exponential behavior exp(— D v /£). As fc^o 
approaches the strong coupling regime within fc M £o < 1, 
however, the rapid oscillation disappears and the expo- 
nential behavior depends on the value of fc M £o- This is 
because the zero energy wave function is spatially ex- 
panded beyond the core region as described in Eq. (|21h>|) 
and implies that the excitation spectrum in the bulk be- 
comes gapless as the topological phase transition point is 
approached. 
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FIG. 2: Lowest eigenenergies as a function of -D v /£o for var- 
ious values of fc M £o i n the case of two- vortex states (N v = 2). 
The symbols are obtained by the fully numerical calculation 
of the BdG equation J5J, while the solid and dashed-dotted 
lines denote exp(—D v /£o)/VD v and exp(— (I — A)D v /£o) with 
A = y/l — (fc M £o) 2 , respectively, which are analytically given 
in Eqs. (|28)l and (J29J) , respectively. 



Then, in order to clarify this splitting and oscillation 
of the Majorana zero modes for all the values of k^ , 
we extend the results derived by Cheng et aZ.J41| to the 
more generic regime beyond the weak coupling limit, i.e., 
kfj,£,o £ (0, oo ) . As two vortices get close to each other with 
decreasing Z3 v /^o ; it is expected that the wave functions 
of the Majorana zero modes bound at each vortex core 
are hybridized with each other, which lifts the degen- 
eracy from zero energy. The variational wave functions 
<p±(p) = [u±(p), v±(p)) T with £ — 0, called the symmetric 
and anti-symmetric states are defined as 



¥>± ; 



V2 



¥?fco,j=i T «¥>fco j=2j 



(22) 



which is valid in the dilute regime of vortices D v 3> £o • 
These variational wave functions fulfill the orthogonal 
condition J <p\<p_ dp = 0. Here, the function <£>fc0j de- 
scribes the wave functions of the zero energy states bound 



at the vortex position Rj=\ 
from Eq. (20|) as 



or R 



■j=2, 



which is obtained 



W,AP) 



e J e 



fe- K +i{Pj)e 



-i(n-l)i 



e t° 



,(23) 



where fij — pj(cos9j,sin8j) denotes the coordinate cen- 
tered at the j-th vortex core, Clj — ^2^j0k{Rj) sums 
up the U(l) phase shift at the j'-th vortex core due to 
the other vortices. Also, the function ft(pj) is described 
in Eq. (J2II) . The particle- hole symmetry obtained from 
Eq. ([7} points out that (fi+(p) with an eigenenergy E is 
associated with tp_{p)=Ti<p\(p) with —E. Note that E 
can have either positive or negative value. Hereafter, we 
refer to the state with u + (v,—) as the symmetric (anti- 
symmetric) state. 



J 



The splitting eigenenergy E+ under the hybridized wave function in Eq. ([22j) obeys the BdG equation, IC(p)<p+(p) = 
E + cp + (p) where K.(p) presented in Eqs. (|A2[) and (|A3|) . corresponding to the low energy approximation of K,(pi,p%) 



in Eq. flU} with T m (k)' 

by 



-k m - With the BdG equations for ipi(p) and ip + (p), the splitting energy E + is now given 



E, 



is V > \(p)^(p)f+(p)dp - | s tp ] + {p)JC{p)ipi{p)dp 
Js^l(p)v+(p)rfp 



(24) 



where £ is defined in the region of x £ [0, oo) and y £ (— oo, oo). The numerator in Eq. (1241) is then explicitly written 
as 

- TpJJ^ J dp [ul(p)\7 2 u 2 (p) + u* 2 {p)\7 2 Ul {p) - c.c] + 4j J dp [u* x {p)Ti.{p)ut{p) + u* 2 (p)n{p)ul{p) - c.c] . 

The first (second) term reduces to the line integral along y at x — by carrying out the integration by parts (by 
employing the Green's theorem). Substituting Eqs. (|2"TT) and (|2"2"|) into Eq. (|M|) . one finds for fc M £o > 1 



E, 



2N 2 b 
M 



dy 



Jo (Vy T +^) Ji (ViF+v 



VtF+v 



£"+ 1 ^y 2 + b 2 



l)j 2 o(V^T¥) 



, Vv 2 +b 2 



(25) 



where J v (z) is the z^-th order Bessel function. The second term in the right hand side of Eq. (|25[) is negligible in the 
weak coupling limit fc^o^^FC^l [4JL|] , while it becomes crucial for the regime near fc M ^o~l- For the strong coupling 
regime within k^o <1, Eq. ((24j) with Eqs. ((22) and (J2TJ reduces to 



E + = 



2N 2 b 



dy 



Vy 2 + b2 



1 



I 



Vv 2 + b2 



- ; > i 



(v^ 2 " 



6 2 



,\/£±E 



(26) 



where ij/(p) is the j/-th order modified Bessel function 
and we introduce 



a± 



'± 



1 



(fc^o) 2 



(27) 



The quantity a± describes the inverse of the wavelength 
of the zero energy state, as seen in Eq. (f2~T|) . In addition, 
the dimensionless variable b is regarded as b = D v a+/2 
in Eq. ([25]) and b = D v a-/2 in Eq. flU}. 

Let us now evaluate the integral in Eq. (|2"5j) for fc^o > 
1. In the leading order of b/£oa+ ^> 1 in Eq. (1231) . i.e., 
the splitting energy is given as 



E 4 



Vo A(l + A 2 ) 
V^f£o (1 + A) 2 
A + cos(Z? v a + ) — A_ sm(D v a^ 
\/D v a + 



' «0 . 



with A± = y / (AyA J TT±A)/(l + A 2 ) and 



A 



h& 



V^(W 



1. 



(28a) 



(28b) 



In the weak coupling limit with k^o « /cf£o 3> 1 and a + ~ 
fcp, this reduces to E + oce~ Dv ^° cos(kpD v + | )/y/kpD v , 
which is consistent with the result in Ref. J4lj . The rapid 
oscillation on the scale of the Fermi wavelength a+ l ^kp 1 
arises from the Bessel function J v (pa+) in the zero en- 
ergy solution described in Eqs. ((20)) and (f2"Tj) . In addi- 
tion, it is found from Eq. (|28|) that regardless of the value 
of fc M £o, the quasiparticle tunneling between neighboring 
vortices splits the degeneracy of the zero energy exponen- 
tially with respect to the ratio of the vortex separation 
D v and the core radius £o ■ This results from the fact that 
the wave function of the Majorana zero modes is local- 
ized within the vortex core ~£o as described in Eqs. (|2"0")) 
and (|2T|) as long as fc^£o > 1. Note that the second term 



of the right-hand side in Eq. (1231) contains the contribu- 
tion e~ Dv '^ without the rapid oscillation. However, the 
leading term turns out to be an order of 0(6~ 3 / 2 ). 

For k^o < 1, the splitting energy is evaluated from 
Eq. ([261) as 



£, 



go A 3 / 2 y/T^A 

y/% fc F £ 



-(1-A) 



«o 



A = a_£ = \/l - (fc/i6 



(29a) 



(29b) 



Here, the deviations from Eq. (f28|) in the regime of 
&m£o > 1 arise from following two points: (i) The dis- 
appearance of the rapid oscillation and (ii) the exponen- 
tial decay factor e~( 1 ~ x * >Dv ^° . Both can be linked to the 
change of the wave functions of the Majorana zero modes 
from J v to I v in Eq. (1211) . In particular, the deviation (ii) 
originates from that the asymptotic form of the modified 
Bessel function I v {z) for a large argument | z\ ^> 1 involves 
the exponential divergence on z and the wave function is 
extended beyond the vortex core region £o- This means 
that for fc M £o > 1, the quasiparticles which occupy the 
Majorana zero modes are able to tunnel between neigh- 
boring vortices over the distance much larger than the 
coherence length £q- In Eq. ([2T))) , the exponential decay 
factor in the resulting expression of E + depends on the 
dimensionless parameter k^o which represents the effec- 
tive pairing interaction. 

Let us now return to Fig. [2] in which the full numeri- 
cal results are compared with the analytical expressions 
derived in Eqs. (|2"g)l and ([2U)h The overall behavior of 
the splitting energy can be well fitted with e~ Dv '^° for 
fc M £o>l and e -(i-A)Ar/£o for fc M £ < 1 given in Eqs. ^ 
and (|2"9l . 

To quantify the self-Hermitian condition of the Majo- 
rana zero modes described in Eq. (|19a[) . u n (p) =u*(p), 
the following quantity is introduced as the relative norm 
between u n (p) and v n (p), 



K^ / [\u n (p)\ 2 -\v n (p)\ 2 ]dp. 



(30) 



This quantity 1Z describes the relative contribution of 
the particle and hole components to the lowest energy 
quasiparticle. Hence, the quantity 1Z is expected to be 
zero for the Majorana zero modes realized in the dilute 
limit of vortices D v >■ £. This is plotted in Fig. [3] for 
several values of fc^£o- This implies that the deviation 
from \u n \ — \v n \ breaks the Majorana nature associated 
with the self-Hermitian relation ry =r\. 



C. Weak coupling BCS regime 

In this subsection, we shall now look more carefully 
into the splitting of the Majorana zero modes in the weak 
coupling regime with fc // £o>l- Here, k^ = k-p is assumed 
for simplicity. The strong coupling regime will be dis- 
cussed in the following subsection. 
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FIG. 3: Relative norm TL between u n (p) and v n (p) for the 
lowest eigenenergy states. All the points correspond to the 
data displayed in Fig. [5] 



Figure Ufa) shows the cross section of the wave func- 
tions \u n (x,y)\ at y = 0, which is obtained from the 
fully numerical calculation of the BdG equation ([HI) for 
kfj,£,o = 6-67 in two- vortex systems N v = 2. The lowest 
eigenenergy at D v = 10.2£ = 34/c F 1 (11.7£o = 39fc F : ) is 
£„/X> o = -6.4xl0~ 6 (-2.9 xl0~ 6 ). Here, it is seen from 
Fig. @|a) that the wave function is exponentially local- 
ized in the range of £o = 6.67fcp centered at the vortex 
cores x = ±D v /2. In addition, the wave function yields 
the rapid oscillation on the scale of the Fermi wavelength 
« 2tt/A:f. 

The wave function displayed in Fig. |U[a) is in good 
agreement with the symmetric <p+(p) in Eq. (|22p which 
results from the hybridization of <pt,j=i and (pg t j=£. Hence, 
the wave function \u n (x, y = 0)| in the case of -D v /£o = 
10.2 of Fig. 0Ja) i s smoothly connected at x = 0, i.e., the 
symmetric solution of the wave function u n (p). In the 
case of -D v /£o = 11.7, however, the wave function of the 
negative energy eigenstate has a node at x = 0, implying 
the anti-symmetric solution of the wave function u n (p). 
The rapid oscillation of the quasiparticle wave function is 
found to affect the energy difference between symmetric 
and anti-symmetric states. 

Hence, it is found that the hybridized wave functions of 
the Majorana zero modes are characterized by two length 
scales £o and Ac 1 — k^ 1 in the weak coupling regime, as 
also seen in Eq. (|2"3"|) . The former length scale correspond- 
ing to the coherence length determines the range of the 
localization of \u(r)\ and \v(r)\ around the core. In the 
same sense as ordinary double well systems 68] , the over- 
lap of the wave function bound at neighboring vortices 
exponentially splits the lowest eigenenergies from zero to 
«± exp(— D v /£ ), as the vortex distance gets close to £ - 
The shorter length scale k^ 1 = kp 1 <c£o in the weak cou- 
pling regime produces the rapid oscillation of the eigenen- 
ergies between the symmetric and anti-symmetric solu- 




FIG. 4: (Color online) (a) Wave functions \u n \ of the lowest 
energy states in fc M £o — 6.67 and N v — 2. Here, the vortex sep- 
aration is set to be -D v /£o = 11.7 (solid line) and 10.2 (dashed- 
dotted line) and the arrows denote the position of the vortices. 
Here, both the solid and dashed-dotted lines correspond to the 
energy states with E/T>o = — 2.9xl0~ 6 and — 6.4xl0~ 6 , respec- 
tively, (b) The lowest eigenenergies are plotted as a function 
of D v /£o £ [10, 18], where the filled (open) circles denote the 
energies of the symmetric (anti-symmetric) state. The inset 
shows the positive eigenenergies in the range of D v /£o € [0, 18] 
with the logarithmic scale. The dashed line in the inset de- 
picts exp(— -D v /£o). 



tions. As we have discussed in Sec. Ill A and described 
in Eq. (f28|) . these length scales are reflected to the en- 
ergy splitting and oscillation of the Majorana zero modes. 
This fact is demonstrated in Fig. 0Jb), where the lowest 
eigenenergies are plotted as a function of -D v /£o- The 
oscillation period is found to be ~ 27rfc F ~ 1 . This agrees 
with the expression of the splitting energy in Eq. (|2"51) for 
the weak coupling regime and the dilute limit of vortices 
D v /£^> 1. The inset in Fig. 0Jb), which shows the same 
data with the logarithmic scale in -D v /£o g [0, 18], also 
demonstrates the exponential dependence of the split- 
ting energy against the vortex separation. Note that the 
splitting of the Majorana zero modes, which is critical 
to the decoherence in the topological quantum computa- 
tion, is also observed in other systems, such as the non- 
abelian quas iholes of the v = | fractional quantum Hall 
state [44|,|45|, Kitaev's honeycomb lattice model (46|, and 
the generic anyon model [47J. 

Figure [3] focuses on the small D v region in the situa- 
tion same as Fig. E^b). In the regime of D v « £qj two 



vortex cores merge and results in one giant hole near the 
center p = 0. Then, it is found from Fig. [5] that the 
rapid oscillation of the eigenenergies with ss k^ disap- 
pears when D v /£ < 2. The corresponding wave function 
\u n {x,y)\ yields a ring shape surrounding the giant vor- 
tex core, as seen in the inset of Fig. El This is contrast 
to the wave function for -D v /£o 3> 1) which consists of 
the double peaks localized at each vortex core. Specif- 
ically, it is seen from the inset of Fig. [5] that \u n {x, y)\ 
with E n /T> = +0.046 at D v /£ = 1.2 includes six phase 
singularities inside the ring at \p\ « £ - This positive 
energy state is continuously connected to the Caroli-de 
Gennes-Matricon core-bound state with azimuthal quan- 
tum number I = 6 at the limit of -D v /£o = [13| . where 
the pair potential results in an axisymmetric vortex with 
winding number 2. Note that the value of £ depends on 
the parameters fc M £o- 



D. Strong coupling regime 

Now, let us turn to the strong coupling regime, fc^o < 
1. Here, the chemical potential remains positive, i.e., 
the Majorana zero mode still survives in the dilute limit 
D v >• £ as discussed in Sec. Ill A. Note that p = or 
fc/j£o = gives rise to the TPT and the further strong 
coupling regime with p, < involves an isotropic trivial 
excitation gap even in the presence of vortices [12|, LL3| • 

The numerical results on the energy splitting are dis- 
played in Fig. [6ja) , where the parameters are set to be 
kfj,£,o — 0.63 and 0.77. Here the numerical results are 
compared with the wave functions <p+(p) based on the 
analytical solution in Eqs. (|21b|) and (|22|) . It is seen that 
the rapid oscillation of the wave function on the scale of 
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FIG. 5: (Color online) Lowest eigenenergies of the core- 
localized states in the range of -D v /£o £ [0, 10]. The param- 
eters are the same as Fig. HJb). The filled and open circles 
denote the symmetric and anti-symmetric states. Three den- 
sity maps describe the quasiparticle wave function |w n (a;, j/)| 
(x,y£ [— 3£o, 3£o]) of the eigenstates classified as the symmet- 
ric state at D v /^o — 1.2, 2.1, and 5.4, respectively. 



Numerical 
Analytical 




W> -40 -30 -20 -10 10 20 30 40 

(a) ~ x/i 

1.5 x 10~ = 

1.0 



e 



0.5, 

-0.51 
-1.0 
-1.5 



*„& = 0.63 

o 



0.77 



O 



symmetric 
O 



y b ■°|2 9 8«§ « ■-■■■• ■■■-■« 



antisymmetric 



10 



15 



(b) 



20 25 



30 



35 



FIG. 6: (Color online) (a) Cross section of the wave function 
\u n (p)\ of the lowest energy states in fc M £o = 0.63 and 0.77 at 
-Dv/£o = 30. The solid lines obtained by the full numerical 
diagonalization of the BdG equation (|14[) are compared with 
the variational wave function based on the analytical solu- 
tion described in Eqs. 1)221) and (|23p (dashed line). Here, the 
vortex separation is fixed to be D v /Co = 30. The solid lines 
correspond to the energy states with E/T>o = — 1.33 x 10 
for fc M £o = 0.63 and E/V = 1.095 x 10~ 7 for fe M f = 0.77, re- 
spectively. The lowest eigenenergies are plotted in (b) as a 
function of D v /£o £ [0,20], where the filled (open) symbols 
denote the energies of the symmetric (anti-symmetric) state 
of u n {p) for ^£0 = 0.63 and 0.77. 



the Fermi wavelength disappears and the wave function 
is extended beyond the core radius £o- This results from 
that the factor e~ Pj ^° in Eq. (|2"Tj) is canceled out by 
the exponential divergence of the modified Bessel func- 
tion for a large p, which depends on the new length scale 
£o/\/l — (kfiS,o) 2 - It is seen from Fig. |6fa) that the re- 
sulting wave functions obtained from the numerical cal- 
culation are found to coincide with the analytic solution 
in Eq. (J2TJ with Eq. (J22|) . 

These modifications of the wave functions arising from 
the strong coupling effect are reflected by the quasipar- 
ticle spectrum displayed in Fig. [6^b). Here, the splitting 
of low-lying eigenenergies is shown as a function of D v in 
the strong coupling regime of fc M £rj = 0.63 and 0.77. The 
negative energy state turns out to be the anti-symmetric 
state of the wave function u(p), while the positive en- 
ergy state is symmetric. In addition, the splitting of the 
eigenenergies is describable with a simple exponential fac- 
tor with the length scale £o/\A — (^»£o) 2 instead of £o- 
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All the results are in good agreement with the expression dispersion, 
analytically derived in Eq. ([%§]) . 



Eq=[Q 



iVv 



k F R' 



(31) 



IV. THREE- VORTEX SYSTEMS 



A. Low-lying edge states in plural vortex systems 

In chiral p-wave supcrfluid with a single vortex, the en- 
ergy spectra of the BdG equation (JSJ consist of two low- 
energy bound states, such as the edge- and core-bound 
states, in addition to the continuum states. Due to the 
odd parity of the pair potential, A(r,— fc) = — A(r,fc), 
the particle incoming to the rigid wall at \r\ = R feels 
the 7r-phase shifted pair potential relative to the outgo- 
ing one. In the weak coupling limit with fc M £ ^ 1) this 
is describable with the one-dimensional Dirac equation 
with mass domain wall [33] , l67j | , which leads to the low 
energy Andreev resonant state bound at domain wall or 
the edge of systems. The dispersion is proportional to 
the azimuthal quantum number £, Eg oc I — (k — l)/2 
when the Cooper pair is in the k x — ik y channel and the 
axisymmetric vortex with the winding number k is as- 
sumed [33] ] . 

Here, we start by demonstrating that the dispersion 
of the edge states in axisymmetric systems is applica- 
ble to the non-axisymmetric situation with plural vor- 
tices. Figure [Tfa) shows the wave function |u„(r)| of the 
lowest edge-bound states in the k x — ik y pairing state 
with N v = 2 vortices. In the weak coupling regime with 
&m?o = 6.67, the wave function |u(p)| is localized within 
the coherence length £,o/R « 0.05 in the vicinity of the 
surface of the system and oscillates with the Fermi wave- 
length. In the strong coupling regime with fc^£o < 1, ( e -9^ 
see the dashed line in Fig. Efa)), the wave function is no 
longer localized within £o and the rapid oscillation on 
the scale of the Fermi wavelength disappears. Further- 
more, it is found from Fig. EJa) and the inset that the 
wave function of the edge state is axially symmetric and 
the phase winding number is well defined along the edge 
of the system, resulting in the axisymmetric expression 
u(p) ~ UQ(p)e l Q e . Hence, the edge-bound states can be 
identified by counting the phase winding number Q of 
the quasiparticle wave function u n (p) or v n {p). 

We display in Figs. [TJb) and [T^c) the energy spectra 
of the edge bound states in the case of N v = 2 and 3 as a 
function of Q. It is here seen that the analytical disper- 
sion relation Egocl — (« + l)/2 is still valid for the non- 
axisymmetric situation with plural vortices. The vortex 
winding number k and the quantum number £ in the dis- 
persion are now replaced to the total number of vortices 
N v and quasi-quantum number QsZ, respectively. Here, 
N v determines the phase winding of the pair potential 
along the closed path near the edge, i.e., A(p, k)(xe lNv6 
at |p|~i2. As shown in Figs. Efb) and[7Jc), the replaced 



is in good agreement with the numerical results of the 
edge states, which allows us to classify the low-lying BdG 
eigenstates with the phase winding Q. It is obvious that 
the lowest energy of Eq can be zero if the total vortex 
number N v is odd. It should be emphasized that this out- 
come is independent of the details of the boundary con- 
dition and the internal information of the system, such as 
D v /£,o- This was demonstrated in our previous work[12j 
that an alternative boundary condition due to a trap po- 
tential which is used to confine atomic clouds does not 
alter the low energy dispersion of the edge bound states 
in vortex- free states. Hence, it is expected that the low- 
est edge bound state for an odd number of N v always 
survives at the zero energy unless the system approaches 
the dense regime of vortices. 
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FIG. 7: (Color online) (a) Wave functions u(x, y = 0) of the 
lowest edge bound state in fc M £o = 6.67 (solid line) and 0.77 
(dashed line) with iV v = 2, where the shape of u and v is 
independent of D v . The inset shows the phase profile of u n (p) 
in the range of x, y £ [-R, R]. The spectrum of the low- lying 
edge state is displayed in (b) and (c) as a function of Q: 
N v = 2 (b) and -/V v =3 (c). The various values of fc^fo are 
plotted here. The solid lines in (b) and (c) correspond to the 
dispersion in Eq. (pTTj) . 
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B. Splitting of Majorana zero modes 

As described above, the edge bound state for three- 
vortex systems contribute to Majorana zero modes as 
well as three core-bound states. Hence, there exist four 
degenerate zero energy states in total, when each vor- 
tex is well separated from each other and the edge. In 
contrast, the low energy spectrum for D v <C£o coincides 
with that in the single "giant" vortex with N v — 1 and 
the winding number n = 3, where the "index theorem" 
ensures the topological protection of one pair of the zero 
energy state [13, bill [32| . Here, we clarify the spectral 
stability of the Majorana zero modes in the intermedi- 
ate regime that the vortex separation and the radius of 
systems are finite but not zero. 

The finiteness of vortex separation and the radius of 
systems lifts the degeneracy of Majorana zero modes 
into four nondegenerate states, where two of them have 
positive energies and the others are found to have the 
corresponding negative energies. This is because of the 
quasiparticle tunneling between the cores and edge. Fig- 
ure[5Ja) shows the wave functions |u„(p)| of the two low- 
est energy states in the weak coupling regime k^o = 4 
along the path A^r B — > C — > A ^ D, numerically ob- 
tained from the BdG equation (fT4)) . The wave function 
in the higher energy state with E n /T>o = 4.6xl0 -7 has 
three peaks at each vortex cores labeled as A, B, and C, 
while the lower energy state with E n /T>o = 2.2 x 10~ 12 
consists of all the contributions from the three core- and 
edge-bound states. It is found that the localization of the 
wave function at the core region is still describable with 
the analytical solution in Eq. (|2"TT) . where they oscillate 
rapidly on the scale of the Fermi wavelength. 

As we have seen in two-vortex systems, the oscillation 
of the wave function gives rise to the oscillation of the 
splitting eigenenergies as a function of the vortex sepa- 
ration. This consequence can be extended to the case 
of three-vortex systems. We present in Fig. [5Jb) the 
low-lying eigenenergies for the three- vortex system with 
fc M £o = 4 as a function of the vortex separation D v /£q. It 
is seen that two of them oscillate on the scale of the Fermi 
wavelength as well as the results in the case of N v — 2. In 
addition, these oscillating eigenenergies are found from 
the inset to obey the exponential splitting with respect 
to Dy/^o, i.e., exp(— D v /£q). The splitting and oscillat- 
ing eigenenergy states consist of three core-bound wave 
functions, as seen in the upper panel of Fig. EJa). The 
state that the edge-bound state contributes always has 
the lower energy than core-bound states and survives in 
the vicinity of the zero energy regardless of D v . It is also 
found from the inset that the lowest eigenenergies are 
embedded around \E n \/T> <10~ 10 , while they gradually 
increase with the rapid oscillation as -D v /£o increases. 
This is due to the quasiparticle tunneling between the 
vortex core and edge as a result of the finite size effect 
with the fixed Rkp = 150. 

Finally, we summarize in Fig.[TJ]the shift of the eigenen- 
ergies as a function of -D v /£o for three-vortex systems. 




FIG. 8: (Color online) (a) Wave functions |u n (p)| of the low- 
est energy states in the case of N v = 3 with £;f£o = 4 and 
D v /£o — 15. As depicted in the upper panel, the horizontal 
axis corresponds to the path along A ~ > B — > C — > A — > D, 
where A, B, and C denote the vortex positions and D is the 
boundary of the system. The upper and bottom panels in (a) 
are referred as the core bound state and edge-core bound state, 
whose energies are E n /T>o =4.6xl0 -7 and 2.2xl0 -12 , respec- 
tively, (b) The lowest eigenenergies are plotted as a function 
of -Dv/Co £ [9, 15], where the filled (open) circles denote the 
energies of the core bound (core-edge bound) state. The inset 
shows the positive eigenenergies in the range of D v /Co £ [0, 18] 
with the logarithmic scale. The dashed line in the inset de- 
picts exp(— D v /£o)- 



For all the data, the upper branches are found to consist 
of the core-bound states without the contribution of the 
edge state as displayed in the upper panel of Fig. EJa). 
The whole behavior of the splitting energy is determined 
by the quasiparticle tunneling between the neighboring 
vortices, which is the same situation as that in the case 
of two- vortex systems. In fact, the upper branches are 
fitted by the exponential factor e~ Dv ^° for k^o > 1 and 
e -\A-(fc^o) 2 £>v/«o for k ^ Q < i as well as Fig.|2]for Nv = 2, 
which implies that the analytical results in Eqs. 
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FIG. 9: (Color online) Lowest eigenenergies as a function of 
-D v /£o for various values of k^t^o in the case of three- vortex 
states (iVy = 3) . The symbols denote the fully numerical solu- 
tion of the BdG equation ([8| and the solid and dashed lines 
depict the analytical results for N v = 2 described in Eqs. (|28[) 
and (1291). 



and (|29[) are applicable to the splitting of Majorana zero 
modes due to the inter-vortex tunneling in three- vortex 
systems. 

The lower branches for each fc M £o in Fig. [9] are identified 
as the core-edge bound states as displayed in the lower 
panel of Fig. Uta). Since the distance between A and 
D becomes short as D v /£q increases, the quasiparticle 
tunneling between the core and edge gives rise to the 
increase of the eigenenergy of the core-edge bound state 
and deviates from the zero energy exponentially. Note 
that this finite-size effect may be negligible as long as the 
system size is macroscopic and the vortices are dilute. 



V. CONCLUDING REMARKS 



scale of kfx v/l — (fc^o) -2 appears in addition to the ex- 
ponential splitting of the eigenenergies. The exponential 
behavior turns out to be characterized uniquely by the ra- 
tio of the vortex separation D v and the coherence length 
of the pair potential £o as e~ Dv ^°. 

In contrast, the rapid oscillation disappears in the 
strong coupling regime close to the topological phase 
transition point at /c M £o = 0. Here, it has been demon- 
strated that the wave function of the zero energy quasi- 
particle bound at a vortex c ore can be ex tended to neigh- 
boring vortices within £ / 



1 — (A; M £o) 2 , which enhances 
the splitting of the Majorana zero energy. These facts 
are reflected by the energy splitting obtained in Eq. ([29]) , 
where the splitting for fc^£o < 1 is uniquely determined by 
the single parameter ^/l — (fc^£o) 2 -Dv/£o- All these be- 
haviors for fc M £ > 1 and fc M £ < 1 have been confirmed by 
numerical diagonalization of the BdG equation (fl4"|) with 
the huge size of the matrix and explained by the drastic 
change of the wave function of the Majorana zero modes 
from the Bessel function to the modified Bessel function. 
Hence, it is concluded that the concrete realization of the 
non-abelian statistics associated with the Majorana zero 
modes requires the neighboring vortices to be separated 
from each other over the length scale £o for fc M £o > 1 and 

£o/v/l-(fc^o) 2 forfc^o<l- 

We have also expanded these arguments into three- 
vortex systems, where one of four zero energy states in di- 
lute limit -D v /£o S> 1 is contributed from the edge-bound 
state. The four degenerate ground states can be differ- 
entiated by the quasiparticle tunneling between vortices 
and/or the edge of the system. The quasiparticle tunnel- 
ing between neighboring vortices gives rise to the split- 
ting and quantum oscillation of two of four degenerate 
zero energy states in the same sense as the case of two- 
vortex systems. In contrast to two- vortex systems, how- 
ever, the other zero energy states composed of the core- 
and edge-bound states is insensitive to the vortex-vortex 
separation but affected by the vortex-edge tunneling, i.e., 
the finite size effect of the system. 



Here, we have investigated the splitting and quantum 
oscillation of Majorana zero modes due to the quasipar- 
ticle tunneling between neighboring vortices in spinless 
chiral p-wave superfluids, based on the analytical and 
numerical calculations of the BdG equation ([M)) . The 
equation which we use in this article contains the non- 
local p-wavc pair potential described in Eq. (flSl) . which 
guarantees all the eigenvalues to be real. In addition, 
the BdG equation reduced to the low energy contains 
the zero energy solution bound at a vortex core when the 
vorticity of the pair potential is odd. 

In two- vortex systems, the analytical expression on 
the splitting eigenenergy of the Majorana zero modes 
has been derived by using the generic solution for the 
zero energy core-bound states. The resulting expressions 
in Eqs. (|28|) and (|29|) tell us that in the weak coupling 
regime within fc^o > 1 the quantum oscillation on the 
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Appendix A: Low-energy approximation on Eq. (|8|) 
and complex eigenvalues 

In this Appendix, we start with the BdG equation (|SJ) 
and the non-local pair potential Eq. © with Eq. (|10l) . 
The expression of the symmetry factor T m (fc) in Eq. (|11[) 
reduces to T m (k) = -r-k m when the low-energy limit is 
taken as k <Cfcp. Then, following the procedure described 
in Refs. 13, 69], the pair potential is expressed with the 
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local function A m (f) and the spatial derivatives P±\ 
=F(9x ± idy) and P = dz as 



A(n,r 2 ) « i^A^WPiMn -r a ) 



(Al) 



Here, we replace i»4 m — > A m - Hence, the matrix in the 
BdG equation (|8]) reduces to the local form as K.(ri,r 2 )~ 
JC(ri)S(ri — r 2 ) with 



JC(r) = 



H (r) n(r) 
-II*(r) -ff *(r) 



H(r) 



1 
2^ 



E 

m=0,±l 



M m (r),P m }. 



(A2a) 



(A2b) 



The resulting BdG equation for the quasiparticles under 
pair potentials in m orbital channels is then given as 



H (r) 

-n*(r) 



n(r) 



u v {r) 
v v (r) 



= E V 



u v {r) 
v v (r) 



(A3) 



The resulting equation turns out to be the eigenvalue 
problem with non-Hcrmitian matrix fC(r), in contrast to 
the original matrix K,(ri,r%) that is Hcrmitian. This is 
in consequence of approximation on the non-local pair 
potential described in Eq. (jAll) , which no longer exhibits 
the p-wave symmetry in Eq. Q . 

The BdG equation (|A3[) simplified to the local form 
is useful for the derivation of the analytic solution of 
the zero energy eigenstates, for instance, as described 
in Eq. (|2"Tj) . Actually, the analytic solutions derived from 
Eq. (|A3|) is in good agreement with the results obtained 
from the numerical diagonalization of the non-local BdG 
equation (|T4")) . Nevertheless, we should notice that the 
simplified equation (IA3|) is not suitable for the investiga- 
tion on accuracy of the zero energy eigenvalues. 

In practice, it is easy to see that the BdG equation (IA3|) 
may contain complex eigenvalues. Assuming the eigen- 
state that yields u v {r) = v*(r), the BdG equation for 
£„eC is rewritten as H {r)u u {r) + H(r)v v (r) = E v u v {r) 
and H (r)uv(r) + U(r)v u (r) = — E*u v (r), leading to 
E v =—E*. Hence, a pure imaginary values E^eC may 
be an eigenvalue as well as the zero energy solution with 
E v = 0, so that they are not distinguishable within the 
numerical diagonalization. Note that the complex eigen- 
values always appear as a pair of E v and — E* because 
of the particle-hole symmetry of the eigenstates derived 
from Eq. (0. 

In Fig. [TU1 we present the lowest eigenvalues in the 
case of three-vortex systems with fc M £o = 6.67 as an ex- 
ample, where the eigenenergies obtained from Eq. (|A3[) 
are compared with those based on the non-local BdG 
equation (1141 with Eq. (|T5|) . It is seen that the real part 
of E v yields abrupt jump at certain value of £> v /£o =8, 
which is the ratio of the vortex separation and coherence 
length. All the eigenvalues for D v /£ < 8 are real, while 
the non- vanishing imaginary part appears in the situa- 
tion beyond D v /£ — 8. In contrast, the imaginary part 
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FIG. 10: (Color online) (a) Real part of the low-lying eigen- 
values E n for ku£o = 6.67 and N v — 3 as a function of the 
vortex separation D v . The filled circles are the same data 
displayed in Fig. IHb) with circles, which result from the di- 
agonalization of Eq. (|14[) with non-local A(n,T"2) in Eq. (|15|l. 
The open squares are the numerical results obtained from 
Eq. (|A3|l . The imaginary part of the corresponding energy 
E n is displayed in (b). 



and abrupt jump of E v never appear in the eigenvalues 
obtained from the BdG equation (fT4| with the non-local 
A(pi,p 2 )- 



Appendix B: Discrete variable representation 

Here, we describe the details about how to numerically 
solve the BdG equation (jHJ with the non-local pair po- 
tential obtained in Eq. (|15[). In order to map the BdG 
equation ((5J into the eigenvalue equation, we apply the 
the discrete variable representation (DVR) [54-56] . The 
details are as follows: First, it is convenient to replace 
the continuous variable p= (x, y) to the N discrete grids 
{xj}j=i,...,jv and {yj}j=i,-,N, where Xj,yj E [-R,R]. 
Here, to interpolate an arbitrary function on these grid 
points, we introduce a set of the N Lagrange polynomials 
within the range g£ [-R, R], 



N 



q-q k 



fM = U^ 



&.i 



Qj - Qk 



which satisfies the conditions, 



fjilk) — Sj,k, 



fj{<l)fk{q)dq 



^j$j,k, 



(Bl) 



(B2a) 



(B2b) 



where q = x,y. The former condition is ensured by the 
Lagrange interpolation with Eq. (IB1|) . The latter can 
be satisfied by setting qj to be the grids obtained by 
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the Gauss-Lobatto quadrature rule [541-156] . where Xj is 
the corresponding weights on qj. The quadrature rule 
turns the integral over the continuous variable q to the 
summation on qj, J G(q)dq « z\,jG{qj)Xj, which gives 
the exact result if the function G(q) is a polynomial of 
degree <2JV-1. 

With those polynomials, let us introduce a set of the 
N orthonormal functions Xj(o) — fi(o)/\\j satisfying 
j-RXj{ ( i)Xk{<l)dq = Sj t k- Then, the eigenfunction in the 
BdG equation ([5]) in the x-y plane are expanded to 



u v {r) 
v v {r) 



N 

E 






(") 



fa 



Xi(x)xj(y)- 



(B3) 



By substituting Eq. (IB3|) to Eq. ((5J) and using the or- 
thonormal condition for Xji the BdG equation ([SJ re- 
duces to the eigenvalue problem with 2iV 2 -dimensional 
Hermite matrix 



H 

At 



A 
-H* 






u v 



(B4) 



O) 



where t/„ and V u are the iV 2 -dimensional vectors of U q 

and Vi? , respectively. H_ and A are the N 2 x iV 2 matrix 
and the elements are obtained as 



{H)ij,i>j> = / dpx t {x)x 3 (y) H o( r )xi> (x)xf (y) 

+ [Vtrap(nj) - M] S i,i'$3,j'' ( B5a ) 



2Af 



ij,i'j' f ij,i'f 



(A)ij,i>j> = dp dp'x l (x)xj{y)^(r,r')x l '(x')xj'{y') 
rj y/X i X j Xi>XfA(p ij , pvji), (B5b) 

where py = (xi,yj). Note that A = —A. In addition, 
the gradient term Tfj V j, is given as |56f 



rpx 



d 2 

dpxiWXj (y) q^Xi' (x)xj' (y) 



hi' E A ' 



dXi(x) dxi'{x) 



dx dx 



(B6a) 



with 



dfM 



dq 



^n 



Qk - q-n 



x=xu & ~ Ik . It - Qn 



■ , for i 7^ k 



■^T-(6i,N — Si,i), iori = k. (B6b) 



Due to the non-locality of the pair potential A(r"i,r2), 
the 2A^ 2 x 2iV 2 matrix elements in Eq. (|B4I) becomes 
dense. The resulting huge and dense matrix is numer- 
ically diagonalized with the shift-invert Lanczos algo- 
rithm [63], which reduces the eigenvalue problem to the 
2iV 2 -dimensional linear equation. Then, the linear equa- 
tion is iteratively solved by using the Krylov subspace 
method, such as the generalized minimal residual method 
with a preconditioner [64|. 
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